﻿using System;
using System.Collections;
using System.Collections.Generic;

namespace Mclib
{
	/// <summary>
	/// This class implements the mapping from a data value to the quantile index in the range (0, 1) and visa versa,
	/// thereby implementing a copula or equalized histogram.
	/// </summary>
	/// <typeparam name="T">The type of the underlying data (usually "double").  Data type must implement IComparable (i.e. it must be sortable).</typeparam>
	public class QuantileList<T> where T : IComparable<T>
	{
		/// <summary>
		/// The list of Sortable values in their original order.
		/// </summary>
		public readonly List<Sortable<T>> OriginalList;
		/// <summary>
		/// The list of Sortable values in their sorted order.
		/// </summary>
		public readonly List<Sortable<T>> SortedList;
		/// <summary>
		/// Creates a QuantileList from a vector.
		/// </summary>
		/// <param name="x">The vector from which to create a quantile list.</param>
		public QuantileList(T[] x)
		{
			int i;
			Sortable<T> s;
			OriginalList = new List<Sortable<T>>(x.Length);
			SortedList = new List<Sortable<T>>(x.Length);
			for (i=0; i<x.Length; i++)
			{
				s = new Sortable<T>(x[i],i);
				OriginalList.Add(s);
				SortedList.Add(s);
			}
			processConstruction();
		}
		/// <summary>
		/// Creates a QuantileList from a matrix.
		/// </summary>
		/// <param name="x">The array from which to create a quantile list from the entries in a single column.</param>
		/// <param name="iCol">The 0-based index of the column of x[,_].</param>
		public QuantileList(T[,] x,int iCol)
		{
			int i;
			int len=x.GetLength(0);
			Sortable<T> s;
			OriginalList = new List<Sortable<T>>(len);
			SortedList = new List<Sortable<T>>(len);
			for (i=0; i<len; i++)
			{
				s = new Sortable<T>(x[i,iCol],i);
				OriginalList.Add(s);
				SortedList.Add(s);
			}
			processConstruction();
		}
		/// <summary>
		/// Used by constructors.
		/// </summary>
		private void processConstruction()
		{
			Sortable<T> s;
			int i,len=SortedList.Count;
			SortedList.Sort();
			SortedList[0].J = 0;
			SortedList[0].Jlow=0;
			for (i=1; i<len; i++)
			{
				s = SortedList[i];
				//	Update J
				s.J = i;
				//	Update Jlow
				if (SortedList[i-1].Value.CompareTo(s.Value)==0)
					s.Jlow = SortedList[i-1].Jlow;
			}
			// Update Jhigh
			SortedList[len-1].Jhigh=len-1;
			for (i=len-1; i>=0; i--)
			{
				s = SortedList[i];
				if (SortedList[i+1].Value.CompareTo(s.Value)==0)
					s.Jhigh=SortedList[i+1].Jhigh;
			}
		}
		/// <summary>
		/// Gets the quantile of a value by locating the nearest values in SortedList via binary search and then
		/// calculating a linear interpolant of the nearest indices.  If val is lower than the first element of SortedList or
		/// greater than the last element, the quantile is extrapolated.
		/// The quantile is always greater than zero (unless val is -Inf) are less than 1 (unless val is +Inf).
		/// </summary>
		/// <param name="val">The value.</param>
		/// <returns>The quantile in the range (0, 1).</returns>
		public float QuantileOf(T val)
		{
			int iSearch = SortedList.BinarySearch(new Sortable<T>(val,-1));
			Sortable<T> s,ss;
			if (iSearch<0)
			{
				//  Exact match was not found.
				//	Get index of the next larger item in the list (bitwise complement).
				iSearch = ~iSearch;
				if (iSearch==0)
				{
					//	Value is out of range (lower)
					ss = SortedList[1];
					s = SortedList[0];
					if (typeof(T)==typeof(double))
					{
						double t = ((double)(object)s.Value-(double)(object)val)/((double)(object)ss.Value-(double)(object)s.Value);
						if (t<1.0)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return (float)((0.5-t/4.0)/(double)SortedList.Count);
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return (float)((0.25/t)/(double)SortedList.Count);
					}
					if (typeof(T)==typeof(float))
					{
						float t = ((float)(object)s.Value-(float)(object)val)/((float)(object)ss.Value-(float)(object)s.Value);
						if (t<1.0f)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return (0.5f-t/4.0f)/(float)SortedList.Count;
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return (0.25f/t)/(float)SortedList.Count;
					}
					if (typeof(T)==typeof(decimal))
					{
						decimal t = ((decimal)(object)s.Value-(decimal)(object)val)/((decimal)(object)ss.Value-(decimal)(object)s.Value);
						if (t<1.0m)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return (float)((0.5m-t/4.0m)/(decimal)SortedList.Count);
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return (float)((0.25m/t)/(decimal)SortedList.Count);
					}
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
				else if (iSearch>=SortedList.Count)
				{
					//	Value is out of range (upper)
					ss = SortedList[SortedList.Count];
					s = SortedList[SortedList.Count-1];
					if (typeof(T)==typeof(double))
					{
						double t = ((double)(object)val-(double)(object)ss.Value)/((double)(object)ss.Value-(double)(object)s.Value);
						if (t<1.0)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return 1.0f-(float)((0.5-t/4.0)/(double)SortedList.Count);
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return 1.0f-(float)((0.25/t)/(double)SortedList.Count);
					}
					if (typeof(T)==typeof(float))
					{
						float t = ((float)(object)val-(float)(object)ss.Value)/((float)(object)ss.Value-(float)(object)s.Value);
						if (t<1.0f)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return 1.0f-(0.5f-t/4.0f)/(float)SortedList.Count;
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return 1.0f-(0.25f/t)/(float)SortedList.Count;
					}
					if (typeof(T)==typeof(decimal))
					{
						decimal t = ((decimal)(object)val-(decimal)(object)ss.Value)/((decimal)(object)ss.Value-(decimal)(object)s.Value);
						if (t<1.0m)
							//	Cast to the range [0.25f-0.5f]/SortedList.Count.
							return 1.0f-(float)((0.5m-t/4.0m)/(decimal)SortedList.Count);
						else
							//	Cast to the range (0, 0.25f]/SortedList.Count
							return 1.0f-(float)((0.25m/t)/(decimal)SortedList.Count);
					}
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
				else
				{
					//	Value is in range
					ss = SortedList[iSearch];
					s = SortedList[iSearch-1];
					if (typeof(T)==typeof(double))
						return (0.5f+(float)s.J+
									(float)(((double)(object)val-(double)(object)s.Value)/((double)(object)ss.Value-(double)(object)s.Value))
								)   /   (float)SortedList.Count;
					if (typeof(T)==typeof(float))
						return (0.5f+(float)s.J+
									((float)(object)val-(float)(object)s.Value)/((float)(object)ss.Value-(float)(object)s.Value)
								)   /   (float)SortedList.Count;
					if (typeof(T)==typeof(decimal))
						return (0.5f+(float)s.J+
									(float)(((decimal)(object)val-(decimal)(object)s.Value)/((decimal)(object)ss.Value-(decimal)(object)s.Value))
								)   /   (float)SortedList.Count;
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
			}
			else
			{
				// Exact match was found.  Return the midpoint of the indexes.
				s = SortedList[iSearch];
				return ((float)(1+s.Jlow+s.Jhigh))/(float)(2*SortedList.Count);
			}
		}
		/// <summary>
		/// Gets the value of a quantile by a floating-point index (0,1) into the SortedList.  The answer is computed by linear interpolation.
		/// </summary>
		/// <param name="quantile">The quantile in the range (0, 1).</param>
		/// <returns>The value.</returns>
		public T ValueOf(float quantile)
		{
			quantile = quantile * (float)SortedList.Count - 0.5f;
			if (quantile>=0.0f)
			{
				float len = (float)(SortedList.Count-1);
				if (quantile<len)
				{
					//	Quantile is in range.
					int fi = (int)quantile;

					if (typeof(T)==typeof(double))
					{
						double imod = (double)quantile-(double)fi;
						double fv = (double)(object)SortedList[fi].Value;
						double cv = (double)(object)SortedList[fi+1].Value;
						return (T)(object)(fv+(cv-fv)*(1.0-imod));
					}
					if (typeof(T)==typeof(float))
					{
						float imod = quantile-(float)fi;
						float fv = (float)(object)SortedList[fi].Value;
						float cv = (float)(object)SortedList[fi+1].Value;
						return (T)(object)(fv+(cv-fv)*(1.0f-imod));
					}
					if (typeof(T)==typeof(decimal))
					{
						decimal imod = (decimal)quantile-(decimal)fi;
						decimal fv = (decimal)(object)SortedList[fi].Value;
						decimal cv = (decimal)(object)SortedList[fi+1].Value;
						return (T)(object)(fv+(cv-fv)*(1.0m-imod));
					}
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
				else if (quantile==len)
					return SortedList[SortedList.Count-1].Value;
				else
				{
					//	Quantile was extrapolated from out of range (upper)
					throw new NotImplementedException("TO DO");
				}
			}
			else
			{
				//	Quantile was extrapolated from out of range (lower)
				if (quantile<=-0.5f)
				{
					if (typeof(T)==typeof(double))
						return (T)(object)double.NegativeInfinity;
					if (typeof(T)==typeof(float))
						return (T)(object)float.NegativeInfinity;
					if (typeof(T)==typeof(decimal))
						return (T)(object)decimal.MinValue;
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
				else
				{
					if (typeof(T)==typeof(double))
					{
						double cv = (double)(object)SortedList[0].Value;
						double fv = (double)(object)SortedList[1].Value;
						fv = cv-(fv-cv);
						if (quantile<-.25f) // (-0.5f, -0.25f)
						{
						}
						else // [-0.25f, 0)
							return (T)(object)(fv+(cv-fv)*(1.0+4.0*(double)quantile));
					}
					if (typeof(T)==typeof(float))
					{
					}
					if (typeof(T)==typeof(decimal))
					{
					}
					throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
				}
			}
			throw new NotImplementedException("Invalid Type.  QuantileList<"+typeof(T).Name+"> not implemented.");
		}
	}
	/// <summary>
	/// This encapsulates an IComparable data type by associating each sample of data with two indices for use in the QuantileList class.
	/// </summary>
	/// <typeparam name="T"></typeparam>
	public class Sortable<T> : IComparable<Sortable<T>>,IComparable<T> where T : IComparable<T>
	{
		/// <summary>
		/// The Sortable Value.
		/// </summary>
		public readonly T Value;
		/// <summary>
		/// The 0-based index of the Sortable Value BEFORE sorting.
		/// </summary>
		public readonly int I;
		/// <summary>
		/// The 0-based index of the Sortable Value AFTER sorting.
		/// </summary>
		public int J=-1;
		/// <summary>
		/// This is equal to J unless this item is equal to adjacent items in the list, in which case this is the lowest index of all these items.
		/// </summary>
		public int Jlow=-1;
		/// <summary>
		/// This is equal to J unless this item is equal to adjacent items in the list, in which case this is the highest index of all these items.
		/// </summary>
		public int Jhigh=-1;
		/// <summary>
		/// Constructs a Sortable item.
		/// </summary>
		/// <param name="Value"></param>
		/// <param name="I"></param>
		public Sortable(T Value,int I)
		{
			this.Value = Value;
			this.I = I;
		}
		#region IComparable Members
		public int CompareTo(Sortable<T> other)
		{
			return Value.CompareTo(other.Value);
		}
		public int CompareTo(T other)
		{
			return Value.CompareTo(other);
		}
		#endregion
	}
}